#import specimen data with outliers removed
specimens_ed <- data.frame(read.csv("../Data/cleaned data/specimen_data_tidy_no-outliers.csv", header=TRUE))
#import species trait data
traits <- data.frame(read.csv("../Data/cleaned data/trait_data_tidy.csv", header=TRUE))
#import phylogeny
tree <- read.nexus("../Data/cleaned data/sergestid_tree_pruned")
# Prep data ------
#Calculate species means (3 outliers removed)
sp_means <- specimens_ed %>%
mutate_if(is.character, as.factor) %>%
group_by(genus_species) %>%
summarise(eye_av = mean(Eye_Diameter),
length_av = mean(Body_Length),
n = n()) %>%
ungroup()
#Merge with species trait data stored in dataframe 'species'
species <- traits %>%
left_join(sp_means, by = "genus_species") %>%
mutate(Organ = factor(Organ, levels = c("pesta","lensed", "unlensed", "none"))) %>%
mutate(species_text = as.factor(str_replace(genus_species, "_", " "))) # Add labels for species text
#check that names match in dataframe and tree
name.check(phy = tree, data = species, data.names = species$genus_species)
#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
species.comp <- comparative.data(data = species, phy = tree, names.col = "genus_species", vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE)
#check for dropped tips or dropped species
species.comp$dropped$tips #phylogeny
species.comp$dropped$unmatched.rows #dataset
Here, we fit PGLS models of body length ~ migration distance, with separate models for lambda = 0 and lambda = 1.
#fit model for eye diameter ~ body length (lambda = 0)
pgls_bodyrange0 <- pgls(length_av ~ Range,
data = species.comp,
lambda = 0.00001, #set lambda to 0
bounds=list(lambda=c(0.00001,0.00001)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_bodyrange0)
par(mfrow = c(1, 1))
#main effects
anova(pgls_bodyrange0)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 0.00, delta = 1.00, kappa = 1.00
##
## Response: length_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Range 1 242.39 242.39 1.4156 0.2539
## Residuals 14 2397.17 171.23
#parameter estimates
summary(pgls_bodyrange0)
##
## Call:
## pgls(formula = length_av ~ Range, data = species.comp, lambda = 1e-05,
## param.CI = 0.95, bounds = list(lambda = c(1e-05, 1e-05)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.925 -13.451 1.005 8.575 18.806
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 0.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.503636 6.267670 4.8668 0.0002493 ***
## Range 0.012781 0.010742 1.1898 0.2539106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.09 on 14 degrees of freedom
## Multiple R-squared: 0.09183, Adjusted R-squared: 0.02696
## F-statistic: 1.416 on 1 and 14 DF, p-value: 0.2539
When lambda = 0, body size is not correlated with the distance of vertical migration.
#fit model for eye diameter ~ body length (lambda = 1)
pgls_bodyrange1 <- pgls(length_av ~ Range,
data = species.comp,
lambda = 1, #set lambda to 1
bounds=list(lambda=c(1,1)),
param.CI = 0.95)
#check model assumptions
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_bodyrange1)
par(mfrow = c(1, 1))
#main effects
anova(pgls_bodyrange1)
## Analysis of Variance Table
## Sequential SS for pgls: lambda = 1.00, delta = 1.00, kappa = 1.00
##
## Response: length_av
## Df Sum Sq Mean Sq F value Pr(>F)
## Range 1 0.08 0.078 4e-04 0.9835
## Residuals 14 2485.19 177.514
#parameter estimates
summary(pgls_bodyrange1)
##
## Call:
## pgls(formula = length_av ~ Range, data = species.comp, lambda = 1,
## param.CI = 0.95, bounds = list(lambda = c(1, 1)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.8326 -15.1206 -5.4750 -0.4621 9.9129
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [Fix] : 1.000
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.20003037 6.40648190 5.6505 5.986e-05 ***
## Range -0.00018434 0.00876772 -0.0210 0.9835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.32 on 14 degrees of freedom
## Multiple R-squared: 3.157e-05, Adjusted R-squared: -0.07139
## F-statistic: 0.000442 on 1 and 14 DF, p-value: 0.9835
When lambda = 1, body size is not correlated with the distance traversed in vertical migration.
Both models agree that body size and vertical migraiton distance are not significantly correlated in this sample of sergestid shrimps.
Here, we plot body size vs. vertical migration distance.
#Give species unique color/shape combinations, with color groups coordinated to phylogenetic subgroups
#reorder factor levels for figure legend (phylo order)
species$species_text <- factor(species$species_text,
levels = c("Deosergestes corniculum",
"Deosergestes henseni",
"Allosergestes pectinatus",
"Allosergestes sargassi",
"Sergestes atlanticus",
"Neosergestes edwardsii",
"Parasergestes vigilax",
"Parasergestes armatus",
"Eusergestes arcticus",
"Gardinerosergia splendens",
"Robustosergia regalis",
"Robustosergia robusta",
"Phorcosergia grandis",
"Sergia tenuiremis",
"Challengerosergia talismani",
"Challengerosergia hansjacobi"))
#make shape pallette
shapes.sp <- c("Deosergestes corniculum" = 21,
"Deosergestes henseni" = 22,
"Allosergestes pectinatus" = 23,
"Allosergestes sargassi" = 24,
"Sergestes atlanticus" = 25,
"Neosergestes edwardsii"= 23,
"Parasergestes vigilax" = 21,
"Parasergestes armatus" = 22,
"Eusergestes arcticus" = 23,
"Gardinerosergia splendens" = 24,
"Robustosergia regalis" = 25,
"Robustosergia robusta" = 21,
"Phorcosergia grandis" = 22,
"Sergia tenuiremis" = 23,
"Challengerosergia talismani" = 24,
"Challengerosergia hansjacobi" = 25)
#sergia/sergestes green purple pallette
cols.sp <- c(#Sergestes group
"Deosergestes corniculum" = "#512E5F",
"Deosergestes henseni" = "#633974",
"Allosergestes pectinatus" = "#76448A",
"Allosergestes sargassi" = "#884EA0",
"Sergestes atlanticus" = "#9B59B6",
"Neosergestes edwardsii" = "#AF7AC5",
"Parasergestes vigilax" = "#C39BD3",
"Parasergestes armatus" = "#D7BDE2",
"Eusergestes arcticus" = "#EBDEF0",
#Sergia group
"Gardinerosergia splendens" = "#ABEBC6",
"Robustosergia regalis" = "#A9DFBF",
"Robustosergia robusta" = "#52BE80",
"Phorcosergia grandis" = "#27AE60",
"Sergia tenuiremis" = "#1E8449",
"Challengerosergia talismani" = "#196F3D",
"Challengerosergia hansjacobi" = "#145A32")
# Make plot
plot_bodyDVM <- ggplot(species, aes(x = Range, y = Body_length, color = species_text, shape = species_text, fill = species_text)) +
geom_point(size = 2, alpha = 1) +
scale_shape_manual(values = shapes.sp, name = "Species") +
scale_color_manual(values = cols.sp, name = "Species") +
scale_fill_manual(values = cols.sp, name = "Species") +
ylab("Body length (mm)") +
xlab("Distance of DVM (m)") +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.text = element_text(face = "italic"))
ggplotly(plot_bodyDVM)